The role of static stress diffusion in the spatio-temporal organization of aftershocks 
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We investigate the spatial distribution of aftershocks and we find that aftershock linear density 
exhibits a maximum, that depends on the mainshock magnitude, followed by a power law decay. 
The exponent controlling the asymptotic decay and the fractal dimensionality of epicenters clearly 
indicate triggering by static stress. The non monotonic behavior of the linear density and its 
dependence on the mainshock magnitude can be interpreted in terms of diffusion of static stress. 
This is supported by the power law growth with exponent H ~ 0.5 of the average main-aftershock 
distance. Implementing static stress diffusion within a stochastic model for aftershock occurrence 
we are able to reproduce aftershock linear density spatial decay, its dependence on the mainshock 
magnitude and its evolution in time. 
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Large earthquakes give rise to a sudden increase of the 
seismic rate in the surrounding area. Aftershocks are of- 
ten observed where mainshocks have increased the static 
Coulomb stress [H, 0, 01 0] and their rate decays in time in 
agreement with state-rate friction laws 0, Q . Aftershocks 
also occur in regions of reduced static stress Q as well 
as at distances up to thousand kms from the mainshock 
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1 31 ] - Dynamic stress related to the pas- 
sage of shock waves, is the most plausible explanation for 
this remote triggering. Many studies, also supported by 
experiments on laboratory fault gouge systems 1J| , have 
recently proposed dynamic stress as the main mechanism 
responsible for aftershock triggering 

0,13, USUI- The 
distribution p(Ar), where Ar is the epicentral distance 
between each aftershock and its related mainshock, rep- 
resents a useful tool to discriminate between triggering 
by static or dynamic stress 



121- 



In both cases, p(Ar) 
is expected to decay asymptotically as Ar _M , where p, is 
related to the fractal dimensionality D of epicenters via 
the relationship p + D — 1 = a with a — 1 or a = 2 for 
dynamic or static stress triggering, respectively. Felzer 
& Brodsky (FB) 0] studied p(Ar) for small and inter- 
mediate mainshock magnitudes, obtaining a pure power 
law decay with an exponent p ~ 1.4. This result, to- 
gether with the estimate D ~ 1, was interpreted in favor 
of dynamic stress triggering aftershocks. In this paper 
we will show that the distribution p(Ar) exhibits a non- 
monotonic behavior, with a power law tail and a max- 
imum depending on the mainshock magnitude that can 
be attributed to a stress diffusion mechanism. 

In our analysis we use the Shearer et al. relocated 
Southern California Catalogue in the years 1981-2005 
with an average uncertainty on the epicentral localiza- 
tion of 0.03 km. We consider all events with magnitude 
m > 2. Mainshocks are identified with the same cri- 
terion used by FB, i.e mainshocks are events separated 
in time and space from larger earthquakes 0]. After- 
shocks are all subsequent events occurring within a cir- 



cular region of radius 100 km centered at the mainshock 
epicenter. In Fig.l we plot p(Ar) for all aftershocks re- 
lated to a mainshock with magnitude m € [M, M + 1 [ 
for M = 2, 3, 4 and for a typical time window of 30 
min post-mainshock, as considered by FB. We find that 
p(Ar) exhibits a maximum at a value of Ar increasing 
with M, followed by a pure power law decay Ar -19 only 
when M = 4. For M = 2,3, conversely, a plateau is 
observed at large distances, Ar > 10km (M = 2) and 
Ar > 30fcm (M = 3), which is related to uncorrelated 
background events. Indeed, p(Ar) can be written as the 
sum p(Ar) = pas (Ar) + ps(Ar), where pAs(Ar) is the 
aftershock density distribution and p#(Ar) oc Ar D ~ l is 
the contribution of background events. Since the after- 
shock number decreases in time whereas background seis- 
micity has a constant rate, p(Ar) ~ p^(Ar) in temporal 
windows sufficiently distant from the mainshock. More 
precisely, we obtain p^(Ar) in temporal widows distant 
more than td = 70 days from the mainshock. Results, 
plotted as open symbols in Fig.l, do not depend on td 
for larger td- For each M, a flat behaviour is obtained 
for Ar > 1 km, implying D ~ 1, in agreement with FB. 
A more precise measurement gives D — 1.03 ± 0.05. The 
value of ps(Ar) depends on M, since it is proportional to 
the number of mainshocks in each class M. This implies 
that ps(Ar) becomes less relevant for larger M and, in 
particular, does not affect the exponent fi — 1.88 ± 0.05 
obtained for M — 4 from Fig.l. For M — 2, 3, conversely, 
the tail of the distribution must be appropriately fitted 
with p(Ar) = ps(Ar) + J 4Ar _M . For Ar > 1 km, the cor- 
relation coefficient provides results consistent with fi = 2 
and excludes \i = 1.4. Hence, the exponent value fi ~ 1.4 
obtained as best fit in the range [0.2 : 50] km (orange line 
in Fig. lb) does not represent the asymptotic decay of 
p(Ar). Similar behavior is obtained for hypocentral dis- 
tances, with small differences only at lengths comparable 
with location errors. 

In order to extend the analysis to larger temporal win- 



dows post-mainshock we use the criterion proposed in 
ref. [l9( to separate aftershocks from background events. 
Given two events with magnitude mi and m 2 with oc- 
currence times t\ < t 2 and locations fi , r 2 , the expected 
number of events inside a circle of radius Ar = \r\ — 
centered in r*i, over a time window T = t 2 — ti is pro- 
portional to n exp (l,2) = CW- b( - mi -^TAr D . Here D = 
1.03, b is the slope of the Gutenberg-Richter magnitude- 
frequency distribution and C — 2.06 10 -11 sec -1 km~ D is 
the average rate of m > 2 earthquakes in the catalog. For 
a given mainshock (fi,ti) each subsequent earthquake 
(f 2 ,t 2 ) with n exp (l,2) < n t /j < 1, where n t /j is a given 
threshold, is highly unexpected and therefore it is consid- 
ered an aftershock. Aftershock number should decay in 
time according to the Omori law, which fixes the value of 
the threshold n t h, in particular we find n t h = 10 -3 . Dif- 
ferent values of D e [1.1, 1.6] provide similar results. This 
criterion allows to discriminate between aftershocks di- 
rectly triggered by the mainshock (first generation) from 
higher order generations, exclu ding eventual effects due 
to aftershock cascading 2(| 21. 122L l23j]. An event 2 is a 
first generation aftershock of the event 1, if in the time in- 
terval ]t±, t 2 [ no event j with n(j, 2) < n(l, 2) is present. 
All the following results are obtained considering only 
first generation aftershocks. No important difference is 
observed if higher order generation aftershocks are in- 
cluded in the analysis. The study of pj^g^Ar) with this 
aftershock selection criterion (Fig. 2) provides results in 
agreement with the previous analysis, i.e. a power law 
decay with an exponent p ~ 2 for all values of M. Fur- 
thermore, curves for different M collapse on the same 
master curve (inset a of Fig. 2) following the scaling 
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with f3 = 0.42 ± 0.02. This result was obtained in ref. 
[l|| using a different mainshock selection criterion. The 
function F(x) is non monotonic and exhibits power law 
behaviour F(x) ~ with p = 1.94 ± 0.04 at large x. 
The collapse of curves with small M on curves with larger 
M, weakly affected by the background seismicity, vali- 
dates the aftershock selection criterion. Fig. 2 confirms 
p ~ 2 supporting the static stress triggering scenario. 

The non-monotonic behaviour of p(Ar) is commonly 
attributed to the violation of the point-source hypothe- 
This implies that seismic sources have a finite 
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extension whose linear size scales with the earthquake 
magnitude Ls{m) = 0.01 10 5m km [HI]. One then com- 
putes pth(Ar) assuming that aftershocks are distributed 
according to a power law from a point randomly chosen 
on the mainshock fault and defining Ar as the distance 
from the center of the mainshock fault. p t h(Ar) follows 
the experimental p{Ar) in the whole spatial range for 
M = 2 (inset (b) in Fig. 2). For larger M, conversely, 
theoretical curves significantly deviates from the experi- 
mental ones. Indeed, curves for different M collapse on 
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FIG. 1: (Color online) The distribution of distances from the 
mainshock (filled symbols) versus Ar for mainshock magni- 
tude m g [M, M + 1[. Aftershocks are events occurring within 
T — 30mm from the mainshock (678, 864, 494 aftershocks for 
M — 2,3,4 respectively). Open symbols represent ps(Ar). 
For M = 4, the power law fit in the range [1 : 100] km gives 
p = 1.88 ±0.05 (dashed blue line in panel c). Magenta curves 
are obtained by adding the experimental ps(Ar) to the nu- 
merical p(Ar). The orange line in panel (b) is the power law 
x~ 1A obtained by FB in the intermediate range [0.2 : 16]km. 
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FIG. 2: (Color online) The distribution of distances from the 
mainshock for M = 2 (circles), M — 3 (squares), M — 4 
(diamonds) and M — 5 (triangles). Aftershocks are events 
occurring within T = 5h from the mainshock. The aftershock 
(mainshock) number is 1065 (12746) for M = 2, 1800 (3410) 
for M = 3, 1425 (349) for M = 4 and 1454 (52) for M = 5. 
Continuous lines are the result of numerical simulations. In 
the inset (a), collapse of the curve is obtained rescaling Ar by 
10^ M according to Eq.(l), with (3 = 0.42. The continuous ma- 
genta line is the theoretical master curve F(x) and the brown 
dashed line shows the asymptotic decay F(x) ~ x~ 2 . In the 
inset (b), comparison of experimental p(Ar) for M — 2,4 
(symbols) with the theoretical predictions p t h{Ar) (continu- 
ous magenta lines). Dashed lines (red M = 2, blue M = 4) 
are the results of the stochastic model simulations. 



the same pure power law decay at distances Ar > Ls(m), 
where the point source hypothesis holds. This implies 
that, even if theoretical curves exhibit a non-monotonic 
behavior, they do not verify the scaling collapse Eq(l). 

The scaling behavior of p(Ar) can be attributed to 
a diffusion process. To this extent, we implement 
static stress diffusion in a stochastic model for seis- 
mic occurrence based on a dynamical scaling assump- 
tion [25|, \2Q, |27j. Within this framework, for a given 
mainshock of magnitude too and an aftershock of mag- 
nitude to, the magnitude difference Am = mo — to, 
Ar and At are not independent variables. More pre- 
cisely, if time is rescaled by a a generic scaling factor A, 
At — ► AAt, the statistical properties are invariant pro- 
vided that Ar -> X H Ar and Ato -> Ato + (1/6) log A, 
where H is a scaling exponent. The scaling relation 
among At, Ar and Ato implies that, for a given main- 
shock of magnitude too, the conditional probability to 
have a magnitude to aftershock at distance Ar after 
a time At, takes the scaling form P(At, Ar, to, too) = 
At -H Gt / ioK^"-" \ G r (£pr)- Under the only assump- 
tion that Gt(y) and G r (x) are normalizable functions, 
one recovers several features of seismic occurrence as the 
GR law, the generalized Omori law, the scaling behavior 
of the intertime distribution [25| . The distribution p{Ar) 
can be obtained by integrating P(At, Ar, to, too)P(too) 
over At and to. The scaling relation for P(At, Ar, to, too) 
and the GR law P(m ) ~ KT f,m ° then give Eq.(l), with 
F(x) oc J °° du dvu~ 1 v~ H G t (u/v)G r (xv~ H ) and (3 
hi I . Assuming the power law decay G T (x) oc for x 
larger than a cut-off xq, F(x) is a non-monotonic func- 
tion with an asymptotic decay F(x) ~ for x ^> 1. 
We therefore implement in the numerical simulations the 
parameters fitted from experimental data, p. = 2 and 
H = 0.47, obtained from j3 = 0.42 and the typical 
value b = 0.9. In particular, following the procedure 
described in (H, we set G t (y) oc (e 1 /^) - 1 + 71)^ 
with the parameters kt — 12.7ft-, 71 = 0.1, b = 0.9, and 
Gr(x) oc (fc r x)~ M for x > x Q with fi = 2, H = 0.48, 
fe r = 5.1 10~ 6 £;TO 2 /sec and x — 3 10~ 3 fcTO/sec z . We 
find that, for all values of M, numerical curves follow 
the experimental ones (Fig. 2). The scaling (1) is fulfilled 
with the numerical F(x) reproducing the experimental 
master curve (inset (a) of Fig. 2). As a further check, we 
add ps(Ar), obtained in Fig.l, to the numerical distri- 
bution p(Ar). Numerical results (Fig.l) very well agree 
with experimental data over the entire spatial range. 

The agreement between experimental and numerical 
results supports the validity of the scaling relation Ar ~ 
At H with H ~ 0.5, which implies that the evolution 
in time of stress is consistent with a diffusion equation. 
More direct evidence of static stress diffusion can be ob- 
tained by the temporal evolution of the main-aftershock 
spatial distance. In particular we compute -Rmax (At) 
(R(At)), i.e. the maximum (average) distance from a 
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FIG. 3: (Color online) The maximum distance RMAx(At) 
for M = 2, 3, 4, 5 (circles, squares, diamonds, triangles) 
from bottom to top. RMAx(At) grows until AtM (Atm = 
16000sec, 2900sec, 690sec, 190sec for M = 5,4,3,2). The 
dashed line is the power law fit RniAx(At) ~ At H with 
H = 0.54 ±0.05 obtained for M = 5 and At < At M - Contin- 
uous lines are the theoretical prediction for LMAx(At). 



mainshock with to G [M,M + 1[, of aftershocks occur- 
ring in the time window [At, At(l + e)]. For all M, 
RMAx(At) exhibits (Fig. 3) a non monotonic behaviour 
with a maximum at a M-dependent typical AtM- AtM 
can be identified as the time when the percentage of 
events identified as aftershocks becomes smaller than the 
90% of the total number of recorded earthquakes. There- 
fore, for At < AtM, n0 significant bias related to the 
aftershock selection procedure is present. In this tem- 
poral regime, similar results are obtained including in 
the analysis all subsequent earthquakes occurring within 
a radius of 100 kms from the mainshock. Fig. 3 shows 
that, for all values of M, Rmax increases for times 
At < At M - For M = 5 where At M = 16000sec, a 
power law regime Rmax (At) ~ At H clearly detected 
with H = 0.54 ± 0.05. On the other hand, the de- 
cay for At > At m originates from a bias introduced 
by the method for aftershock selection. The condition 
n exp (l, 2) < nth, indeed, implies that aftershocks are only 
events occurring within a given temporal-magnitude re- 
gion and, in particular, all events occurring at distances 
larger than L M Ax(At) oc 10"*/° (At)- 1 / are not con- 
sidered as aftershocks. The tails of Rmax (At) are con- 
sistent with a pure power law decay At~ x / D in agreement 
with the analytical expression for LMAx(At) (Fig. 3). 

Further indication of diffusion can be obtained in the 
regime At > AtM by considering the average distance 
i?(At) inside a region of radius L sup . This can be eval- 
uated as R(At) = f^ 3up dArArp(Ar)/ j^ 3up dArp(Ar), 
using the decay p(Ar) oc (Ar+K )~ 2 obtained from Fig. 2 
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FIG. 4: (Color online) The average spatial distance R(At) 
from the mainshock (black circles) of aftershocks occurring 
in the time window [At, At(l + e)] with e = 0.02 versus At. 
The initial growth R(At) ~ At 0,47 (magenta line) is con- 
sistent with the diffusion behaviour. The decay at larger 
times is related to the upper cut-off L sup . Continuous red 
and dashed green curves are the theoretical prediction Eq.(2): 
Green curves correspond to K — 1.8 io°- 42 ( M - 5 > km fitted 
from Fig. 2. Red curves correspond to K = 0.018 At 0,47 km, 
where At is measured in seconds, obtained from the numerical 
model. 



According to the previous analysis L s 
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At < Atjvf and L sup — L max (At) when At > ialm- 
We introduce in the above equation K = BAt H with 
B = 0.018km/sec H and H = 0.47, obtained fr om the nu- 
merical simulations. Fig. 4 shows that for all M, without 
any further parameter tuning, the theoretical prediction 
(2) reproduces experimental results in the whole time 
range. In Fig. 4 we also plot Eq.2 assuming a constant 
K, obtained as the best fit from Fig. 2. In this case, the 
theoretical R(At) (dashed lines in Fig. 4) overestimates 
the experimental R(At) at small At, whereas it some- 
how underestimates it at larger times. Previous analyses 
20l 21 , 2^, 23[ have obtained a smaller value of the dif- 
fusion exponent, H ~ 0.1. The basic differences with 
our study is that in previous analyses aftershocks have 
not been classified according to the mainshock magnitude 
and distances significantly smaller than the mainshock 
fault length have been included in the analysis. Inter- 
estingly, McKernon and Main [22J recover H ~ 0.5 at 
very large distances, where the point source hypothesis 
is recovered. 

In conclusion, we have shown that static stress is the 
main mechanism responsible for aftershock occurrence. 
Indeed, by properly taking into account background seis- 
micity, p(Ar) exhibits the scaling behavior (1) with the 
power law decay expected within the static stress trig- 
gering scenario. Moreover, the very good agreement of 
the theoretical prediction (2) with the numerical results 
and experimental data indicates that the aftershock spa- 



tial organization evolves in time according to a diffu- 
sion equation. Migration of aftershocks [28| is often ob- 
served and interpreted within different contexts, includ- 
ing state/rate friction 0, Q , viscoelastic relaxation pro - 
H, IS 31 1 and aftershock cascading 20l 21. 22 



cess 

In the present study, the latter mechanism can be dis- 
carded, since only aftershocks directly triggered by the 
mainshock have been considered. The estimated value 
B = 0.018km/ sec H predicts, on average, a post seismic 
stress change over a region of about 10 2 km in 7 years. 
This is consistent with simulations of 3d viscoelastic post 
seismic relaxation after the 1992 Landers earthquake 



[1 

[2 

[3 

[4 
[5 
[6 
[7] 
[8] 
[9] 

[io; 
in 

[12 
[13 

[14 
[15 

[16 

[17 

[18 

[19 

[20 

[21 
[22 

[23 
[24 
[25 

[26 

[27 

[28 
[29 



Reasenberg, P.A., Simpson, R.W., Science 255, 1687 
(1992). 

King, G.C.P., Stein, R.S., Lin, J., Bull. Seism. Soc. Am. 
84, 935 (1994). 

Hardebeck, J.L., Nazareth, J.J., Hauksson, E., J. Geo- 
phys. Res. -Solid Earth 103, 24427 (1998). 
Wyss, M., Wiemer, S., Science 290, 1334 (2000). 
Dieterich, JH. J. Geophys. Res 99, 2601 (1994). 
Stein R.S., Nature 402, pp. 605-609 (1999). 
Parsons, T., J. Geophys. Res. 107 2199,1, (2002). 
Hill DP. et al., Science 260, 1617, (1993). 
Stark,M.A., Davis, S.D., Geophys. Res. Lett. 23, 945, 
(1996). 

Brodsky, E.E., Karakostas, V., Kanamori, H., Geophys. 
Res. Lett. 27, 2741, (2000). 

Gomberg,J., Reasenberg, P.A., Bodin, P., Harris, R.A., 
Nature 411, 462, (2001). 

Eberhart-Phillips, D., et al., Science 300, 1113, (2003). 
Gomberg, J., Bodin, P., Larson, K., Dragert, H., Nature 
427, 621, (2004). 

Johnson, P, Jia, X., Nature 437, 871, (2005). 

Kilb, D., Gomberg, J., Bodin, P., Nature 408, 570, 

(2000). 

Kilb, D., Gomberg, J., Bodin, P., J. Geophys. Res. 107 
doi: 10. 1029 /2001JB0002002 , (2002) . 
Felzer, K.R., Brodsky, E.E., Nature 441, 735, (2006). 
Shearer, P., Hauksson, E., Lin, G., Bull. Seismol. Soc. 
Am. 95, 904, (2005). 

Baiesi, M., Paczuski, M., Phys. Rev. E 69, 0661061 
(2004). 

Helmstetter, A., Sornette, D., Phys. Rev. E , 66 , 061104, 
(2002). 

Hue, M., Main, I.G., J. Geophys. Res., 108, 2324, (2003). 
McKernon, C, Main, I.G., J. Geophys. Res., 110, 
B05S05, (2005). 

Marsan, D. & Lengline, O., Science, 319, 1076, (2008). 
Kagan Y.Y., BSSA 92, 641 (2002). 

Lippiello, E., Godano, C, de Arcangelis, L., Phys. Rev. 
Lett. 98, 098501 (2007). 

Lippiello, E., de Arcangelis, L., Godano, C, Phys. Rev. 
Lett. 100, 038501 (2008). 

Lippiello, E., Bottiglieri, M., Godano, C, de Arcangelis, 

L., Geophys. Res. Lett. 34, L23301 (2007). 

Mogi, K., 1968,Bull. Earth. Res. Inst. Univ. Tokyo 46, 57, 

(1968). 

Rydelek, P.A., Sacks, I.S., Nature, 336, 234, (1988). 



■5 



[30] Freed, A.M., Lin, J., Nature 441, 180, (2001). 
[31] Jonsson, S., Nature Geoscience 1, 136 (2008). 



